Part 1: Importing and Combining Timeseries Datasets¶
In [1]:
#import libraries
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
pd.set_option('display.float_format', lambda x: '%.4f' % x)
import seaborn as sns
sns.set_context('paper', font_scale=1.3)
sns.set_style('white')
import warnings
warnings.filterwarnings('ignore')
from time import time
import matplotlib.ticker as tkr
from scipy import stats
from statsmodels.tsa.stattools import adfuller
from sklearn import preprocessing
from statsmodels.tsa.stattools import pacf
%matplotlib inline
import math
import keras
from keras.models import Sequential
from keras.layers import Dense
from keras.layers import LSTM
from keras.layers import Dropout
from keras.layers import *
from sklearn.preprocessing import MinMaxScaler
from sklearn.metrics import mean_squared_error
from sklearn.metrics import mean_absolute_error
from keras.callbacks import EarlyStopping
In [2]:
#uncomment if using the outdated dataset
# #create a for loop that reads all sheets from the excel file and updates the column headers
# x = [2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021]
# for i in x:
# globals()[f"platts{i}"] = pd.read_excel('Platts Historical.xlsx', sheet_name=f"{i}")
# globals()[f"platts{i}"].columns = globals()[f"platts{i}"].iloc[0]
# globals()[f"platts{i}"] = globals()[f"platts{i}"].iloc[2:]
# globals()[f"platts{i}"].reset_index(inplace=True, drop=True)
In [3]:
# #contenate all of the created dataframes into one big dataframe
# df_comp = pd.concat([platts2009, platts2010, platts2011, platts2012, platts2013, platts2014, platts2015, platts2016, platts2017, platts2018, platts2019, platts2020, platts2021], ignore_index=True, axis=0)
# df_comp.head()
In [4]:
# df_comp = df_comp.rename(columns={np.nan: 'Date'})
# df_comp.head()
In [5]:
# #save the one big dataframe into a csv file for checkpoint and later calling
# df_comp.to_csv('mergedplatts.csv', index = False)
In [6]:
#uncomment if using the outdated dataset
# df = pd.read_csv('mergedplatts.csv')
# df.head()
In [7]:
from google.colab import drive
drive.mount('/content/drive')
In [8]:
df = pd.read_excel("/content/drive/My Drive/Decisions Data Dump/Platts/Petroleum Data.xlsx")
In [9]:
df1 = df.copy()
In [10]:
df.head()
Out[10]:
In [11]:
df['Date']= pd.to_datetime(df['Date'])
df['year']= df['Date'].apply(lambda x: x.year)
df['quarter']= df['Date'].apply(lambda x: x.quarter)
df['month']= df['Date'].apply(lambda x: x.month)
df['day']= df['Date'].apply(lambda x: x.day)
In [12]:
df.head()
Out[12]:
In [13]:
df= df.reset_index(drop=True)
df['weekday']= df.apply(lambda row: row['Date'].weekday(), axis=1)
df['weekday']= (df['weekday']<5).astype(int) #to capture weekdays
print(df.shape)
print(df.Date.min())
print(df.Date.max())
df.tail(5)
Out[13]:
In [14]:
print('The time series starts from: ', df.Date.min())
print('The time series ends on: ', df.Date.max())
Part 2: EDA¶
In [15]:
# gasoline = df[["Date","PGAEY00", "year", "quarter", "month", "day"]].dropna().rename(columns={"PGAEY00": "Price"})
# print(gasoline.shape)
# gasoline.tail()
In [16]:
# diesel = df[["Date", "AAPPF00", "year", "quarter", "month", "day"]].dropna().rename(columns={"AAPPF00": "Price"})
# print(diesel.shape)
# diesel.head()
Plot¶
In [17]:
# g1=gasoline.loc[:,['Date','Price']]
# g1.set_index('Date', inplace=True)
# d1=diesel.loc[:,['Date','Price']]
# d1.set_index('Date', inplace=True)
In [18]:
df.set_index('Date', inplace=True)
In [19]:
df.Gasoline.plot(figsize = (25,5), label='Gasoline', color='#e58a3d', linewidth=1)
df.Diesel.plot(figsize = (25,5), label='Diesel', color='#3f5b9c', linewidth=1)
df['Jet Fuel'].plot(figsize = (25,5), label='Jet Fuel', color='#c24e5d', linewidth=1)
plt.ylabel('Price')
plt.legend()
plt.tight_layout()
plt.title('Diesel vs Gasoline vs Jet Fuel Prices', size=20)
sns.despine(top=True)
plt.show();
Boxplot of yearly vs quarterly Price¶
In [20]:
plt.figure(figsize=(15,10))
plt.subplot(3,2,1)
plt.subplots_adjust(wspace=0.2)
sns.boxplot(x='year', y='Diesel', data=df)
plt.xlabel('year')
plt.title('Box plot of Diesel Yearly Price')
sns.despine(left=True)
plt.tight_layout()
plt.subplot(3,2,2)
sns.boxplot(x='quarter', y='Diesel', data=df)
plt.xlabel('quarter')
plt.title('Box plot of Diesel Quarterly Price')
sns.despine(left=True)
plt.tight_layout()
plt.subplot(3,2,3)
plt.subplots_adjust(wspace=0.2)
sns.boxplot(x='year', y='Gasoline', data=df)
plt.xlabel('year')
plt.title('Box plot of Gasoline Yearly Price')
sns.despine(left=True)
plt.tight_layout()
plt.subplot(3,2,4)
sns.boxplot(x='quarter', y='Gasoline', data=df)
plt.xlabel('quarter')
plt.title('Box plot of Gasoline Quarterly Price')
sns.despine(left=True)
plt.tight_layout();
plt.subplot(3,2,5)
plt.subplots_adjust(wspace=0.2)
sns.boxplot(x='year', y='Jet Fuel', data=df)
plt.xlabel('year')
plt.title('Box plot of Jet Fuel Yearly Price')
sns.despine(left=True)
plt.tight_layout()
plt.subplot(3,2,6)
sns.boxplot(x='quarter', y='Jet Fuel', data=df)
plt.xlabel('quarter')
plt.title('Box plot of Jet Fuel Quarterly Price')
sns.despine(left=True)
plt.tight_layout();
Price distribution (QQ Plot)¶
In [21]:
plt.figure(figsize=(14,20))
plt.subplot(3,2,1)
plt.subplots_adjust(wspace=0.2)
df.Diesel.hist(bins=50, color='#3f5b9c')
plt.title('Diesel Price Distribution')
plt.subplot(3,2,2)
stats.probplot(df.Diesel, plot=plt)
df.Diesel.describe().T
plt.subplot(3,2,3)
plt.subplots_adjust(wspace=0.2)
df.Gasoline.hist(bins=50, color='#e58a3d')
plt.title('Gasoline Price Distribution')
plt.subplot(3,2,4)
stats.probplot(df.Gasoline, plot=plt)
df.Gasoline.describe().T
plt.subplot(3,2,5)
plt.subplots_adjust(wspace=0.2)
df['Jet Fuel'].hist(bins=50, color='#c24e5d')
plt.title('Jet Fuel Price Distribution')
plt.subplot(3,2,6)
stats.probplot(df['Jet Fuel'], plot=plt)
df['Jet Fuel'].describe().T
Out[21]:
Resampling¶
In [22]:
fig= plt.figure(figsize=(18,16))
fig.subplots_adjust(hspace=.4)
ax1= fig.add_subplot(6,1,1)
ax1.plot(df.Diesel.resample('D').mean(), linewidth=2, color='#3f5b9c')
ax1.set_title('Mean Diesel Price resampled over day')
ax2= fig.add_subplot(6,1,2, sharex=ax1)
ax2.plot(df.Diesel.resample('W').mean(), linewidth=2, color='#3f5b9c')
ax2.set_title('Mean Diesel Price resampled over week')
ax2.tick_params(axis='both', which='major')
ax3= fig.add_subplot(6,1,3, sharex=ax1)
ax3.plot(df.Diesel.resample('M').mean(), linewidth=2, color='#3f5b9c')
ax3.set_title('Mean Diesel Price resampled over month')
ax3.tick_params(axis='both', which='major')
ax4= fig.add_subplot(6,1,4, sharex=ax1)
ax4.plot(df.Diesel.resample('Q').mean(), linewidth=2, color='#3f5b9c')
ax4.set_title('Mean Diesel Price resampled over quarter')
ax4.tick_params(axis='both', which='major')
ax5= fig.add_subplot(6,1,5, sharex=ax1)
ax5.plot(df.Diesel.resample('A').mean(), linewidth=2, color='#3f5b9c')
ax5.set_title('Mean Diesel Price resampled over year')
ax5.tick_params(axis='both', which='major')
ax6= fig.add_subplot(6,1,6, sharex=ax1)
ax6.plot(df.Diesel.resample('SM').mean(), linewidth=2, color='#3f5b9c')
ax6.set_title('Mean Diesel Price resampled over semi-month end')
ax6.tick_params(axis='both', which='major')
In [23]:
fig= plt.figure(figsize=(18,16))
fig.subplots_adjust(hspace=.4)
ax1= fig.add_subplot(6,1,1)
ax1.plot(df.Gasoline.resample('D').mean(), linewidth=2, color='#e58a3d')
ax1.set_title('Mean Gasoline Price resampled over day')
ax2= fig.add_subplot(6,1,2, sharex=ax1)
ax2.plot(df.Gasoline.resample('W').mean(), linewidth=2, color='#e58a3d')
ax2.set_title('Mean Gasoline Price resampled over week')
ax2.tick_params(axis='both', which='major')
ax3= fig.add_subplot(6,1,3, sharex=ax1)
ax3.plot(df.Gasoline.resample('M').mean(), linewidth=2, color='#e58a3d')
ax3.set_title('Mean Gasoline Price resampled over month')
ax3.tick_params(axis='both', which='major')
ax4= fig.add_subplot(6,1,4, sharex=ax1)
ax4.plot(df.Gasoline.resample('Q').mean(), linewidth=2, color='#e58a3d')
ax4.set_title('Mean Gasoline Price resampled over quarter')
ax4.tick_params(axis='both', which='major')
ax5= fig.add_subplot(6,1,5, sharex=ax1)
ax5.plot(df.Gasoline.resample('A').mean(), linewidth=2, color='#e58a3d')
ax5.set_title('Mean Gasoline Price resampled over year')
ax5.tick_params(axis='both', which='major')
ax6= fig.add_subplot(6,1,6, sharex=ax1)
ax6.plot(df.Gasoline.resample('SM').mean(), linewidth=2, color='#e58a3d')
ax6.set_title('Mean Gasoline Price resampled over semi-month end')
ax6.tick_params(axis='both', which='major')
In [24]:
fig= plt.figure(figsize=(18,16))
fig.subplots_adjust(hspace=.4)
ax1= fig.add_subplot(6,1,1)
ax1.plot(df['Jet Fuel'].resample('D').mean(), linewidth=2, color='#c24e5d')
ax1.set_title('Mean Jet Fuel Price resampled over day')
ax2= fig.add_subplot(6,1,2, sharex=ax1)
ax2.plot(df['Jet Fuel'].resample('W').mean(), linewidth=2, color='#c24e5d')
ax2.set_title('Mean Jet Fuel Price resampled over week')
ax2.tick_params(axis='both', which='major')
ax3= fig.add_subplot(6,1,3, sharex=ax1)
ax3.plot(df['Jet Fuel'].resample('M').mean(), linewidth=2, color='#c24e5d')
ax3.set_title('Mean Jet Fuel Price resampled over month')
ax3.tick_params(axis='both', which='major')
ax4= fig.add_subplot(6,1,4, sharex=ax1)
ax4.plot(df['Jet Fuel'].resample('Q').mean(), linewidth=2, color='#c24e5d')
ax4.set_title('Mean Jet Fuel Price resampled over quarter')
ax4.tick_params(axis='both', which='major')
ax5= fig.add_subplot(6,1,5, sharex=ax1)
ax5.plot(df['Jet Fuel'].resample('A').mean(), linewidth=2, color='#c24e5d')
ax5.set_title('Mean Jet Fuel Price resampled over year')
ax5.tick_params(axis='both', which='major')
ax6= fig.add_subplot(6,1,6, sharex=ax1)
ax6.plot(df['Jet Fuel'].resample('SM').mean(), linewidth=2, color='#c24e5d')
ax6.set_title('Mean Jet Fuel Price resampled over semi-month end')
ax6.tick_params(axis='both', which='major')
Plot mean Price grouped by year, quarter, month, and day¶
In [25]:
plt.figure(figsize=(14,8))
plt.subplot(2,2,1)
df.groupby('year').Diesel.agg('mean').plot(color='#3f5b9c', linewidth=3)
plt.xlabel('')
plt.title('Mean Diesel Price by Year')
plt.subplot(2,2,2)
df.groupby('quarter').Diesel.agg('mean').plot(color='#3f5b9c', linewidth=3)
plt.xlabel('')
plt.title('Mean Diesel Price by Quarter')
plt.subplot(2,2,3)
df.groupby('month').Diesel.agg('mean').plot(color='#3f5b9c', linewidth=3)
plt.xlabel('')
plt.title('Mean Diesel Price by Month')
plt.subplot(2,2,4)
df.groupby('day').Diesel.agg('mean').plot(color='#3f5b9c', linewidth=3)
plt.xlabel('')
plt.title('Mean Diesel Price by Day');
In [26]:
plt.figure(figsize=(14,8))
plt.subplot(2,2,1)
df.groupby('year').Gasoline.agg('mean').plot(color='#e58a3d', linewidth=3)
plt.xlabel('')
plt.title('Mean Gasoline Price by Year')
plt.subplot(2,2,2)
df.groupby('quarter').Gasoline.agg('mean').plot(color='#e58a3d', linewidth=3)
plt.xlabel('')
plt.title('Mean Gasoline Price by Quarter')
plt.subplot(2,2,3)
df.groupby('month').Gasoline.agg('mean').plot(color='#e58a3d', linewidth=3)
plt.xlabel('')
plt.title('Mean Gasoline Price by Month')
plt.subplot(2,2,4)
df.groupby('day').Gasoline.agg('mean').plot(color='#e58a3d', linewidth=3)
plt.xlabel('')
plt.title('Mean Gasoline Price by Day');
In [27]:
plt.figure(figsize=(14,8))
plt.subplot(2,2,1)
df.groupby('year')['Jet Fuel'].agg('mean').plot(color='#c24e5d', linewidth=3)
plt.xlabel('')
plt.title('Mean Jet Fuel Price by Year')
plt.subplot(2,2,2)
df.groupby('quarter')['Jet Fuel'].agg('mean').plot(color='#c24e5d', linewidth=3)
plt.xlabel('')
plt.title('Mean Jet Fuel Price by Quarter')
plt.subplot(2,2,3)
df.groupby('month')['Jet Fuel'].agg('mean').plot(color='#c24e5d', linewidth=3)
plt.xlabel('')
plt.title('Mean Jet Fuel Price by Month')
plt.subplot(2,2,4)
df.groupby('day')['Jet Fuel'].agg('mean').plot(color='#c24e5d', linewidth=3)
plt.xlabel('')
plt.title('Mean Jet Fuel Price by Day');
Part 3: Statistical Normality and Stationarity Tests¶
Normality¶
In [28]:
stat, p = stats.normaltest(df.Diesel)
print('Statistics = %.3f, p = %.3f' % (stat,p))
alpha = 0.05
if p > alpha:
print('Data looks Gaussian (fail to reject H0)')
else:
print('Data does not look Gaussian (reject H0)')
In [29]:
print(df.Diesel.describe())
In [30]:
sns.distplot(df.Diesel, color='#3f5b9c');
print('Kurtosis: {}'.format(stats.kurtosis(df.Diesel)))
print('Skewness: {}'.format(stats.skew(df.Diesel)))
In [31]:
stat, p = stats.normaltest(df.Gasoline)
print('Statistics = %.3f, p = %.3f' % (stat,p))
alpha = 0.05
if p > alpha:
print('Data looks Gaussian (fail to reject H0)')
else:
print('Data does not look Gaussian (reject H0)')
In [32]:
print(df.Gasoline.describe())
In [33]:
sns.distplot(df.Gasoline, color='#e58a3d');
print('Kurtosis: {}'.format(stats.kurtosis(df.Gasoline)))
print('Skewness: {}'.format(stats.skew(df.Gasoline)))
In [34]:
stat, p = stats.normaltest(df['Jet Fuel'])
print('Statistics = %.3f, p = %.3f' % (stat,p))
alpha = 0.05
if p > alpha:
print('Data looks Gaussian (fail to reject H0)')
else:
print('Data does not look Gaussian (reject H0)')
In [35]:
print(df['Jet Fuel'].describe())
In [36]:
sns.distplot(df['Jet Fuel'], color='#c24e5d');
print('Kurtosis: {}'.format(stats.kurtosis(df['Jet Fuel'])))
print('Skewness: {}'.format(stats.skew(df['Jet Fuel'])))
Dickey-Fuller Test¶
In [37]:
df1['Date']= pd.to_datetime(df1['Date'])
In [38]:
df1.head()
Out[38]:
In [39]:
def test_stationarity(timeseries):
rolmean = timeseries.rolling(window=30).mean()
rolstd = timeseries.rolling(window=30).std()
plt.figure(figsize=(14,5))
sns.despine(left=True)
orig = plt.plot(timeseries, color='blue',label='Original')
mean = plt.plot(rolmean, color='red', label='Rolling Mean')
std = plt.plot(rolstd, color='black', label = 'Rolling Std')
plt.legend(loc='best'); plt.title('Rolling Mean & Standard Deviation')
plt.show()
print ('<Results of Dickey-Fuller Test>')
dftest = adfuller(timeseries, autolag='AIC')
dfoutput = pd.Series(dftest[0:4],
index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
for key,value in dftest[4].items():
dfoutput['Critical Value (%s)'%key] = value
print(dfoutput)
In [40]:
test_stationarity(df1['Diesel'])
In [41]:
test_stationarity(df1.Gasoline)
In [42]:
test_stationarity(df1['Jet Fuel'])
Part 4: Converting Non Stationary Time Series to Stationary using Differencing¶
In [43]:
diesel = df1[['Date','Diesel']]
gasoline = df1[['Date','Gasoline']]
jet = df1[['Date','Jet Fuel']]
In [44]:
#transform to stationary
diesel["diff_1"] = diesel["Diesel"].diff(periods=1)
diesel.head()
Out[44]:
In [45]:
test_stationarity(diesel["diff_1"].dropna())
In [46]:
#transform to stationary
gasoline["diff_1"] = gasoline["Gasoline"].diff(periods=1)
gasoline.head()
Out[46]:
In [47]:
test_stationarity(gasoline["diff_1"].dropna())
In [48]:
#transform to stationary
jet["diff_1"] = jet["Jet Fuel"].diff(periods=1)
jet.head()
Out[48]:
In [49]:
test_stationarity(jet["diff_1"].dropna())
Part 5: LSTM¶
Gasoline LSTM¶
In [50]:
gasoline.set_index('Date', inplace=True)
gasoline.head()
Out[50]:
In [51]:
gasoline1=gasoline.dropna()
print(gasoline1.shape)
gasoline1.head()
Out[51]:
In [52]:
#convert column to array
np.reshape(gasoline1.diff_1.values, (-1,1))
Out[52]:
In [53]:
dataset= gasoline1.diff_1.values #numpy.ndarray
dataset= dataset.astype('float32')
dataset= np.reshape(dataset,(-1,1))
In [54]:
scaler= MinMaxScaler(feature_range=(0,1))
dataset= scaler.fit_transform(dataset)
dataset #compressed
Out[54]:
In [55]:
len(dataset)
Out[55]:
In [56]:
# train_size = int(len(dataset) * 0.80)
# test_size = len(dataset) - train_size
train_size=2905
train, test = dataset[0:train_size,:], dataset[train_size:len(dataset),:]
In [57]:
print(train.shape, test.shape)
In [58]:
#convert an array of values into a dataset matrix
def create_dataset(dataset, look_back=1):
#creates window pockkets
X,Y = [], []
for i in range(len(dataset)-look_back-1):
a=dataset[i:(i+look_back),0]
X.append(a)
Y.append(dataset[i + look_back,0])
return np.array(X), np.array(Y)
In [59]:
#reshape into X=t and Y=t+1
look_back = 30
X_train, Y_train = create_dataset(train, look_back)
X_test, Y_test= create_dataset(test, look_back)
In [60]:
print(X_train.shape, Y_train.shape, X_test.shape, Y_test.shape)
In [61]:
X_test
Out[61]:
In [62]:
#reshape input to be [samples, time steps, features]
X_train= np.reshape(X_train, (X_train.shape[0], 1, X_train.shape[1]))
X_test = np.reshape(X_test, (X_test.shape[0],1, X_test.shape[1]))
In [63]:
X_test
Out[63]:
In [64]:
print(X_train.shape, Y_train.shape, X_test.shape, Y_test.shape)
In [65]:
#creating layers
model = Sequential()
model.add(LSTM(200, input_shape=(X_train.shape[1], X_train.shape[2]))) #LSTM layer with 100 memory units
model.add(Dropout(0.2))
model.add(Dense(1)) # no activation applied
model.compile(loss='mean_squared_error', optimizer='adam')
history = model.fit(X_train, Y_train, epochs=20, batch_size=70, validation_data=(X_test, Y_test),
callbacks=[EarlyStopping(monitor='val_loss', patience=20)], verbose=1, shuffle=False)
# Training Phase
model.summary()
In [66]:
train_predict= model.predict(X_train)
test_predict= model.predict(X_test)
#invert predictions
train_predict= scaler.inverse_transform(train_predict)
Y_train=scaler.inverse_transform([Y_train])
test_predict= scaler.inverse_transform(test_predict)
Y_test= scaler.inverse_transform([Y_test])
print('Train Mean Absolute Error:', mean_absolute_error(Y_train[0], train_predict[:,0]))
print('Train Root Mean Squared Error:',np.sqrt(mean_squared_error(Y_train[0], train_predict[:,0])))
print('Test Mean Absolute Error:', mean_absolute_error(Y_test[0], test_predict[:,0]))
print('Test Root Mean Squared Error:',np.sqrt(mean_squared_error(Y_test[0], test_predict[:,0])))
In [67]:
print(train_predict.shape, test_predict.shape)
In [68]:
plt.figure(figsize=(8,4))
plt.plot(history.history['loss'], label='Train Loss')
plt.plot(history.history['val_loss'], label='Test Loss')
plt.title('model loss')
plt.ylabel('loss')
plt.xlabel('epochs')
plt.legend(loc='upper right')
plt.show();
In [69]:
aa=[x for x in range(60)]
plt.figure(figsize=(8,4))
plt.plot(aa, Y_test[0][:60], marker='.', label="actual")
plt.plot(aa, test_predict[:,0][:60], 'r', label="prediction")
# plt.tick_params(left=False, labelleft=True) #remove ticks
plt.tight_layout()
sns.despine(top=True)
plt.subplots_adjust(left=0.07)
plt.ylabel('Price', size=15)
plt.xlabel('Time step', size=15)
plt.legend(fontsize=15)
plt.show();
In [70]:
aa=[x for x in range(300)]
plt.figure(figsize=(8,4))
plt.plot(aa, Y_train[0][:300], marker='.', label="actual")
plt.plot(aa, train_predict[:,0][:300], 'r', label="prediction")
# plt.tick_params(left=False, labelleft=True) #remove ticks
plt.tight_layout()
sns.despine(top=True)
plt.subplots_adjust(left=0.07)
plt.ylabel('Price', size=15)
plt.xlabel('Time step', size=15)
plt.legend(fontsize=15)
plt.show();
In [71]:
from sklearn.metrics import r2_score
print('Train R2 Score: ', r2_score(Y_train[0], train_predict[:,0]))
print('Test R2 Score: ', r2_score(Y_test[0], test_predict[:,0]))
In [72]:
# !pip install pmdarima
In [73]:
from pmdarima.utils import diff_inv
def inv_diff (df_orig_column,df_diff_column, periods):
# Generate np.array for the diff_inv function - it includes first n values(n =
# periods) of original data & further diff values of given periods
value = np.array(df_orig_column[:periods].tolist()+df_diff_column[periods:].tolist())
# Generate np.array with inverse diff
inv_diff_vals = diff_inv(value, periods,1 )[periods:]
return inv_diff_vals
In [74]:
gasoline['r_diff'] = inv_diff(gasoline.Gasoline, gasoline.diff_1, 1)
gasoline.head()
Out[74]:
In [75]:
# print(r_train.index.min(), r_train.index.max())
In [76]:
# print(r_train.index.min(), r_train.index.max())
In [77]:
r_test = gasoline.r_diff[train_size+2:-(look_back)]
len(r_test)
Out[77]:
In [78]:
print(r_test.index.min(), r_test.index.max())
In [79]:
#convert test_predict to series
test_predict_s = pd.Series(test_predict.reshape(-1))
In [80]:
test_predict_s
Out[80]:
In [81]:
#use def inv_diff to get inverse-diff of predicted test values
r_predict = inv_diff(df['Gasoline'][train_size:], test_predict_s, 1)
In [82]:
r_predict_s = pd.Series(r_predict)
r_predict_s.index = r_test.index
In [83]:
print(r_test.index.min(), r_test.index.max())
# print(r_train.index.min(), r_train.index.max())
print(r_predict_s.index.min(), r_predict_s.index.max())
In [84]:
df['Gasoline'][:train_size].plot(figsize = (25,5), label='train', color='#e58a3d')
df['Gasoline'][train_size:].plot(figsize = (25,5), label='test', color='#e53d98')
r_predict_s.plot(figsize = (25,5), label='prediction', color='#e53d98', linestyle='--')
plt.title('LSTM Gasoline Results', size=24)
plt.legend()
plt.show()
In [85]:
print('Test MSE after Inverse Diff:', mean_absolute_error(r_test, r_predict_s))
print('Test RMSE after Inverse Diff:',np.sqrt(mean_absolute_error(r_test, r_predict_s)))
print('Test R2 Score after Inverse Diff: ', r2_score(r_test, r_predict_s))
In [86]:
lstm_gas = r_predict_s.copy()
In [87]:
len(lstm_gas)
Out[87]:
Diesel LSTM¶
In [88]:
diesel.set_index('Date', inplace=True)
diesel.head()
Out[88]:
In [89]:
diesel1=diesel.dropna()
print(diesel1.shape)
diesel1.head()
Out[89]:
In [90]:
#convert column to array
np.reshape(diesel1.diff_1.values, (-1,1))
Out[90]:
In [91]:
dataset= diesel1.diff_1.values #numpy.ndarray
dataset= dataset.astype('float32')
dataset= np.reshape(dataset,(-1,1))
In [92]:
scaler= MinMaxScaler(feature_range=(0,1))
dataset= scaler.fit_transform(dataset)
dataset #compressed
Out[92]:
In [93]:
len(dataset)
Out[93]:
In [94]:
# train_size = int(len(dataset) * 0.80)
# test_size = len(dataset) - train_size
train_size = 2905
train, test = dataset[0:train_size,:], dataset[train_size:len(dataset),:]
In [95]:
print(train.shape, test.shape)
In [96]:
#convert an array of values into a dataset matrix
def create_dataset(dataset, look_back=1):
#creates window pockkets
X,Y = [], []
for i in range(len(dataset)-look_back-1):
a=dataset[i:(i+look_back),0]
X.append(a)
Y.append(dataset[i + look_back,0])
return np.array(X), np.array(Y)
In [97]:
#reshape into X=t and Y=t+1
look_back = 30
X_train, Y_train = create_dataset(train, look_back)
X_test, Y_test= create_dataset(test, look_back)
In [98]:
print(X_train.shape, Y_train.shape, X_test.shape, Y_test.shape)
In [99]:
X_test
Out[99]:
In [100]:
#reshape input to be [samples, time steps, features]
X_train= np.reshape(X_train, (X_train.shape[0], 1, X_train.shape[1]))
X_test = np.reshape(X_test, (X_test.shape[0],1, X_test.shape[1]))
In [101]:
X_test
Out[101]:
In [102]:
print(X_train.shape, Y_train.shape, X_test.shape, Y_test.shape)
In [103]:
#creating layers
model = Sequential()
model.add(LSTM(200, input_shape=(X_train.shape[1], X_train.shape[2]))) #LSTM layer with 100 memory units
model.add(Dropout(0.2))
model.add(Dense(1)) # no activation applied
model.compile(loss='mean_squared_error', optimizer='adam')
history = model.fit(X_train, Y_train, epochs=20, batch_size=70, validation_data=(X_test, Y_test),
callbacks=[EarlyStopping(monitor='val_loss', patience=20)], verbose=1, shuffle=False)
# Training Phase
model.summary()
In [104]:
train_predict= model.predict(X_train)
test_predict= model.predict(X_test)
#invert predictions
train_predict= scaler.inverse_transform(train_predict)
Y_train=scaler.inverse_transform([Y_train])
test_predict= scaler.inverse_transform(test_predict)
Y_test= scaler.inverse_transform([Y_test])
print('Train Mean Absolute Error:', mean_absolute_error(Y_train[0], train_predict[:,0]))
print('Train Root Mean Squared Error:',np.sqrt(mean_squared_error(Y_train[0], train_predict[:,0])))
print('Test Mean Absolute Error:', mean_absolute_error(Y_test[0], test_predict[:,0]))
print('Test Root Mean Squared Error:',np.sqrt(mean_squared_error(Y_test[0], test_predict[:,0])))
In [105]:
print(train_predict.shape, test_predict.shape)
In [106]:
plt.figure(figsize=(8,4))
plt.plot(history.history['loss'], label='Train Loss')
plt.plot(history.history['val_loss'], label='Test Loss')
plt.title('model loss')
plt.ylabel('loss')
plt.xlabel('epochs')
plt.legend(loc='upper right')
plt.show();
In [107]:
aa=[x for x in range(60)]
plt.figure(figsize=(8,4))
plt.plot(aa, Y_test[0][:60], marker='.', label="actual")
plt.plot(aa, test_predict[:,0][:60], 'r', label="prediction")
# plt.tick_params(left=False, labelleft=True) #remove ticks
plt.tight_layout()
sns.despine(top=True)
plt.subplots_adjust(left=0.07)
plt.ylabel('Price', size=15)
plt.xlabel('Time step', size=15)
plt.legend(fontsize=15)
plt.show();
In [108]:
aa=[x for x in range(300)]
plt.figure(figsize=(8,4))
plt.plot(aa, Y_train[0][:300], marker='.', label="actual")
plt.plot(aa, train_predict[:,0][:300], 'r', label="prediction")
# plt.tick_params(left=False, labelleft=True) #remove ticks
plt.tight_layout()
sns.despine(top=True)
plt.subplots_adjust(left=0.07)
plt.ylabel('Price', size=15)
plt.xlabel('Time step', size=15)
plt.legend(fontsize=15)
plt.show();
In [109]:
from sklearn.metrics import r2_score
print('Train R2 Score: ', r2_score(Y_train[0], train_predict[:,0]))
print('Test R2 Score: ', r2_score(Y_test[0], test_predict[:,0]))
In [110]:
# !pip install pmdarima
In [111]:
from pmdarima.utils import diff_inv
def inv_diff (df_orig_column,df_diff_column, periods):
# Generate np.array for the diff_inv function - it includes first n values(n =
# periods) of original data & further diff values of given periods
value = np.array(df_orig_column[:periods].tolist()+df_diff_column[periods:].tolist())
# Generate np.array with inverse diff
inv_diff_vals = diff_inv(value, periods,1 )[periods:]
return inv_diff_vals
In [112]:
diesel['r_diff'] = inv_diff(diesel.Diesel, diesel.diff_1, 1)
diesel.head()
Out[112]:
In [113]:
r_train = diesel.r_diff[1:train_size+1]
len(r_train)
Out[113]:
In [114]:
print(r_train.index.min(), r_train.index.max())
In [115]:
print(r_train.index.min(), r_train.index.max())
In [116]:
r_test = diesel.r_diff[train_size+2:-(look_back)]
len(r_test)
Out[116]:
In [117]:
print(r_test.index.min(), r_test.index.max())
In [118]:
#convert test_predict to series
test_predict_s = pd.Series(test_predict.reshape(-1))
In [119]:
test_predict_s
Out[119]:
In [120]:
#use def inv_diff to get inverse-diff of predicted test values
r_predict = inv_diff(r_test, test_predict_s, 1)
In [121]:
r_predict_s = pd.Series(r_predict)
r_predict_s.index = r_test.index
In [122]:
print(r_test.index.min(), r_test.index.max())
print(r_train.index.min(), r_train.index.max())
print(r_predict_s.index.min(), r_predict_s.index.max())
In [123]:
df['Diesel'][:train_size].plot(figsize = (25,5), label='train', color='#3f5b9c')
df['Diesel'][train_size:].plot(figsize = (25,5), label='test', color='#9c3f5b')
r_predict_s.plot(figsize = (25,5), label='prediction', linestyle='--', color='#9c3f5b')
plt.title('LSTM Diesel Results', size=24)
plt.legend()
plt.show()
In [124]:
print('Test MSE after Inverse Diff:', mean_absolute_error(r_test, r_predict_s))
print('Test RMSE after Inverse Diff:',np.sqrt(mean_absolute_error(r_test, r_predict_s)))
print('Test R2 Score after Inverse Diff: ', r2_score(r_test, r_predict_s))
In [125]:
lstm_diesel = r_predict_s.copy()
In [126]:
len(lstm_diesel)
Out[126]:
Jet Fuel LSTM¶
In [127]:
jet.set_index('Date', inplace=True)
jet.head()
Out[127]:
In [128]:
jet1=jet.dropna()
print(jet1.shape)
jet1.head()
Out[128]:
In [129]:
#convert column to array
np.reshape(jet1.diff_1.values, (-1,1))
Out[129]:
In [130]:
dataset= jet1.diff_1.values #numpy.ndarray
dataset= dataset.astype('float32')
dataset= np.reshape(dataset,(-1,1))
In [131]:
scaler= MinMaxScaler(feature_range=(0,1))
dataset= scaler.fit_transform(dataset)
dataset #compressed
Out[131]:
In [132]:
len(dataset)
Out[132]:
In [133]:
# train_size = int(len(dataset) * 0.80)
# test_size = len(dataset) - train_size
train_size = 2905
train, test = dataset[0:train_size,:], dataset[train_size:len(dataset),:]
In [134]:
print(train.shape, test.shape)
In [135]:
#convert an array of values into a dataset matrix
def create_dataset(dataset, look_back=1):
#creates window pockkets
X,Y = [], []
for i in range(len(dataset)-look_back-1):
a=dataset[i:(i+look_back),0]
X.append(a)
Y.append(dataset[i + look_back,0])
return np.array(X), np.array(Y)
In [136]:
#reshape into X=t and Y=t+1
look_back = 30
X_train, Y_train = create_dataset(train, look_back)
X_test, Y_test= create_dataset(test, look_back)
In [137]:
print(X_train.shape, Y_train.shape, X_test.shape, Y_test.shape)
In [138]:
X_test
Out[138]:
In [139]:
#reshape input to be [samples, time steps, features]
X_train= np.reshape(X_train, (X_train.shape[0], 1, X_train.shape[1]))
X_test = np.reshape(X_test, (X_test.shape[0],1, X_test.shape[1]))
In [140]:
X_test
Out[140]:
In [141]:
print(X_train.shape, Y_train.shape, X_test.shape, Y_test.shape)
In [142]:
#creating layers
model = Sequential()
model.add(LSTM(200, input_shape=(X_train.shape[1], X_train.shape[2]))) #LSTM layer with 100 memory units
model.add(Dropout(0.2))
model.add(Dense(1)) # no activation applied
model.compile(loss='mean_squared_error', optimizer='adam')
history = model.fit(X_train, Y_train, epochs=20, batch_size=70, validation_data=(X_test, Y_test),
callbacks=[EarlyStopping(monitor='val_loss', patience=20)], verbose=1, shuffle=False)
# Training Phase
model.summary()
In [143]:
train_predict= model.predict(X_train)
test_predict= model.predict(X_test)
#invert predictions
train_predict= scaler.inverse_transform(train_predict)
Y_train=scaler.inverse_transform([Y_train])
test_predict= scaler.inverse_transform(test_predict)
Y_test= scaler.inverse_transform([Y_test])
print('Train Mean Absolute Error:', mean_absolute_error(Y_train[0], train_predict[:,0]))
print('Train Root Mean Squared Error:',np.sqrt(mean_squared_error(Y_train[0], train_predict[:,0])))
print('Test Mean Absolute Error:', mean_absolute_error(Y_test[0], test_predict[:,0]))
print('Test Root Mean Squared Error:',np.sqrt(mean_squared_error(Y_test[0], test_predict[:,0])))
In [144]:
print(train_predict.shape, test_predict.shape)
In [145]:
plt.figure(figsize=(8,4))
plt.plot(history.history['loss'], label='Train Loss')
plt.plot(history.history['val_loss'], label='Test Loss')
plt.title('model loss')
plt.ylabel('loss')
plt.xlabel('epochs')
plt.legend(loc='upper right')
plt.show();
In [146]:
aa=[x for x in range(60)]
plt.figure(figsize=(8,4))
plt.plot(aa, Y_test[0][:60], marker='.', label="actual")
plt.plot(aa, test_predict[:,0][:60], 'r', label="prediction")
# plt.tick_params(left=False, labelleft=True) #remove ticks
plt.tight_layout()
sns.despine(top=True)
plt.subplots_adjust(left=0.07)
plt.ylabel('Price', size=15)
plt.xlabel('Time step', size=15)
plt.legend(fontsize=15)
plt.show();
In [147]:
aa=[x for x in range(300)]
plt.figure(figsize=(8,4))
plt.plot(aa, Y_train[0][:300], marker='.', label="actual")
plt.plot(aa, train_predict[:,0][:300], 'r', label="prediction")
# plt.tick_params(left=False, labelleft=True) #remove ticks
plt.tight_layout()
sns.despine(top=True)
plt.subplots_adjust(left=0.07)
plt.ylabel('Price', size=15)
plt.xlabel('Time step', size=15)
plt.legend(fontsize=15)
plt.show();
In [148]:
from sklearn.metrics import r2_score
print('Train R2 Score: ', r2_score(Y_train[0], train_predict[:,0]))
print('Test R2 Score: ', r2_score(Y_test[0], test_predict[:,0]))
In [149]:
# !pip install pmdarima
In [150]:
from pmdarima.utils import diff_inv
def inv_diff (df_orig_column,df_diff_column, periods):
# Generate np.array for the diff_inv function - it includes first n values(n =
# periods) of original data & further diff values of given periods
value = np.array(df_orig_column[:periods].tolist()+df_diff_column[periods:].tolist())
# Generate np.array with inverse diff
inv_diff_vals = diff_inv(value, periods,1 )[periods:]
return inv_diff_vals
In [151]:
jet['r_diff'] = inv_diff(jet['Jet Fuel'], jet.diff_1, 1)
jet.head()
Out[151]:
In [152]:
r_train = jet.r_diff[1:train_size+1]
len(r_train)
Out[152]:
In [153]:
print(r_train.index.min(), r_train.index.max())
In [154]:
print(r_train.index.min(), r_train.index.max())
In [155]:
r_test = jet.r_diff[train_size+2:-(look_back)]
len(r_test)
Out[155]:
In [156]:
print(r_test.index.min(), r_test.index.max())
In [157]:
#convert test_predict to series
test_predict_s = pd.Series(test_predict.reshape(-1))
In [158]:
test_predict_s
Out[158]:
In [159]:
#use def inv_diff to get inverse-diff of predicted test values
r_predict = inv_diff(r_test, test_predict_s, 1)
In [160]:
r_predict_s = pd.Series(r_predict)
r_predict_s.index = r_test.index
In [161]:
print(r_test.index.min(), r_test.index.max())
print(r_train.index.min(), r_train.index.max())
print(r_predict_s.index.min(), r_predict_s.index.max())
In [162]:
df['Jet Fuel'][:train_size].plot(figsize = (25,5), label='train', color='#c24e5d')
df['Jet Fuel'][train_size:].plot(figsize = (25,5), label='test', color='#5dc24e')
r_predict_s.plot(figsize = (25,5), label='prediction', linestyle='--', color='#5dc24e')
plt.title('LSTM Jet Fuel Results', size=24)
plt.legend()
plt.show()
In [163]:
print('Test MSE after Inverse Diff:', mean_absolute_error(r_test, r_predict_s))
print('Test RMSE after Inverse Diff:',np.sqrt(mean_absolute_error(r_test, r_predict_s)))
print('Test R2 Score after Inverse Diff: ', r2_score(r_test, r_predict_s))
In [164]:
lstm_jet = r_predict_s.copy()
len(lstm_jet)
Out[164]:
Part 6: Stacked LSTM¶
Gasoline Stacked LSTM¶
In [165]:
# gasoline.set_index('Date', inplace=True)
gasoline.head()
Out[165]:
In [166]:
df1=gasoline.Gasoline
df1
Out[166]:
In [167]:
scaler= MinMaxScaler(feature_range=(0,1))
df1=scaler.fit_transform(np.array(df1).reshape(-1,1))
df1 #compressed
Out[167]:
In [168]:
len(df1)
Out[168]:
In [169]:
# train_size = int(len(df1) * 0.80)
# test_size = len(df1) - train_size
train_size = 2906
train, test = df1[0:train_size,:], df1[train_size:len(df1),:]
In [170]:
print(train.shape, test.shape)
In [171]:
#convert an array of values into a dataset matrix
def create_dataset(dataset, look_back=1):
#creates window pockkets
X,Y = [], []
for i in range(len(dataset)-look_back-1):
a=dataset[i:(i+look_back),0]
X.append(a)
Y.append(dataset[i + look_back,0])
return np.array(X), np.array(Y)
In [172]:
#reshape into X=t and Y=t+1
look_back = 30
X_train, Y_train = create_dataset(train, look_back)
X_test, Y_test= create_dataset(test, look_back)
In [173]:
print(X_train.shape, Y_train.shape, X_test.shape, Y_test.shape)
In [174]:
# reshape input to be [samples, time steps, features] which is required for LSTM
X_train =X_train.reshape(X_train.shape[0],X_train.shape[1] , 1)
X_test = X_test.reshape(X_test.shape[0],X_test.shape[1] , 1)
In [175]:
print(X_train.shape, Y_train.shape, X_test.shape, Y_test.shape)
In [176]:
model=Sequential()
model.add(LSTM(100,return_sequences=True,input_shape=(30,1)))
model.add(LSTM(100,return_sequences=True))
model.add(LSTM(100))
model.add(Dense(1))
model.compile(loss='mean_squared_error',optimizer='adam')
In [177]:
model.summary()
In [178]:
history = model.fit(X_train, Y_train, epochs=20, batch_size=70, validation_data=(X_test, Y_test),
callbacks=[EarlyStopping(monitor='val_loss', patience=20)], verbose=1, shuffle=False)
In [179]:
train_predict= model.predict(X_train)
test_predict= model.predict(X_test)
#invert predictions
train_predict= scaler.inverse_transform(train_predict)
Y_train=scaler.inverse_transform([Y_train])
test_predict= scaler.inverse_transform(test_predict)
Y_test= scaler.inverse_transform([Y_test])
print('Train Mean Absolute Error:', mean_absolute_error(Y_train[0], train_predict[:,0]))
print('Train Root Mean Squared Error:',np.sqrt(mean_squared_error(Y_train[0], train_predict[:,0])))
print('Test Mean Absolute Error:', mean_absolute_error(Y_test[0], test_predict[:,0]))
print('Test Root Mean Squared Error:',np.sqrt(mean_squared_error(Y_test[0], test_predict[:,0])))
In [180]:
plt.figure(figsize=(8,4))
plt.plot(history.history['loss'], label='Train Loss')
plt.plot(history.history['val_loss'], label='Test Loss')
plt.title('model loss')
plt.ylabel('loss')
plt.xlabel('epochs')
plt.legend(loc='upper right')
plt.show();
In [181]:
aa=[x for x in range(60)]
plt.figure(figsize=(8,4))
plt.plot(aa, Y_test[0][:60], marker='.', label="actual")
plt.plot(aa, test_predict[:,0][:60], 'r', label="prediction")
# plt.tick_params(left=False, labelleft=True) #remove ticks
plt.tight_layout()
sns.despine(top=True)
plt.subplots_adjust(left=0.07)
plt.ylabel('Price', size=15)
plt.xlabel('Time step', size=15)
plt.legend(fontsize=15)
plt.show();
In [182]:
aa=[x for x in range(300)]
plt.figure(figsize=(8,4))
plt.plot(aa, Y_train[0][:300], marker='.', label="actual")
plt.plot(aa, train_predict[:,0][:300], 'r', label="prediction")
# plt.tick_params(left=False, labelleft=True) #remove ticks
plt.tight_layout()
sns.despine(top=True)
plt.subplots_adjust(left=0.07)
plt.ylabel('Price', size=15)
plt.xlabel('Time step', size=15)
plt.legend(fontsize=15)
plt.show();
In [183]:
from sklearn.metrics import r2_score
print('Train R2 Score: ', r2_score(Y_train[0], train_predict[:,0]))
print('Test R2 Score: ', r2_score(Y_test[0], test_predict[:,0]))
In [184]:
len(test)
Out[184]:
In [185]:
x_input=test[70:].reshape(1,-1)
x_input.shape
Out[185]:
In [186]:
temp_input=list(x_input)
temp_input=temp_input[0].tolist()
In [187]:
# demonstrate prediction for next 10 days
from numpy import array
lst_output=[]
n_steps=30
i=0
while(i<30):
if(len(temp_input)>30):
#print(temp_input)
x_input=np.array(temp_input[1:])
print("{} day input {}".format(i,x_input))
x_input=x_input.reshape(1,-1)
x_input = x_input.reshape((1, n_steps, 1))
#print(x_input)
yhat = model.predict(x_input, verbose=0)
print("{} day output {}".format(i,yhat))
temp_input.extend(yhat[0].tolist())
temp_input=temp_input[1:]
#print(temp_input)
lst_output.extend(yhat.tolist())
i=i+1
else:
x_input = x_input.reshape((1, n_steps,1))
yhat = model.predict(x_input, verbose=0)
print(yhat[0])
temp_input.extend(yhat[0].tolist())
print(len(temp_input))
lst_output.extend(yhat.tolist())
i=i+1
print(lst_output)
In [188]:
pred = pd.Series( (v[0] for v in lst_output) )
In [189]:
day_new=np.arange(1,31)
day_pred=np.arange(31,61)
In [190]:
len(df1)
Out[190]:
In [191]:
plt.plot(day_new,scaler.inverse_transform(df1[2976:]))
plt.plot(day_pred,scaler.inverse_transform(lst_output))
Out[191]:
In [192]:
print(Y_test.shape, Y_train.shape, test_predict.shape)
In [193]:
gasoline.index[:train_size-31]
Out[193]:
In [194]:
test_predict.shape
Out[194]:
In [195]:
plt.figure(figsize=(25,5))
plt.plot(df['Gasoline'][:train_size], label='train', color='#e58a3d')
plt.plot(df['Gasoline'][train_size:], label='test', color='#e53d98')
plt.plot(gasoline.index[-69:], test_predict[:,0], label='prediction', color='#e53d98', linestyle='--')
plt.title('Stacked LSTM Gasoline Results', size=24)
plt.legend()
plt.show()
In [196]:
stlstm_gas = test_predict[:,0].copy()
In [197]:
len(stlstm_gas)
Out[197]:
Diesel Stacked LSTM¶
In [198]:
diesel.head()
Out[198]:
In [199]:
df1=diesel.Diesel
df1
Out[199]:
In [200]:
scaler= MinMaxScaler(feature_range=(0,1))
df1=scaler.fit_transform(np.array(df1).reshape(-1,1))
df1 #compressed
Out[200]:
In [201]:
len(df1)
Out[201]:
In [202]:
# train_size = int(len(df1) * 0.80)
# test_size = len(df1) - train_size
train_size = 2906
train, test = df1[0:train_size,:], df1[train_size:len(df1),:]
In [203]:
print(train.shape, test.shape)
In [204]:
#convert an array of values into a dataset matrix
def create_dataset(dataset, look_back=1):
#creates window pockkets
X,Y = [], []
for i in range(len(dataset)-look_back-1):
a=dataset[i:(i+look_back),0]
X.append(a)
Y.append(dataset[i + look_back,0])
return np.array(X), np.array(Y)
In [205]:
#reshape into X=t and Y=t+1
look_back = 30
X_train, Y_train = create_dataset(train, look_back)
X_test, Y_test= create_dataset(test, look_back)
In [206]:
print(X_train.shape, Y_train.shape, X_test.shape, Y_test.shape)
In [207]:
# reshape input to be [samples, time steps, features] which is required for LSTM
X_train =X_train.reshape(X_train.shape[0],X_train.shape[1] , 1)
X_test = X_test.reshape(X_test.shape[0],X_test.shape[1] , 1)
In [208]:
print(X_train.shape, Y_train.shape, X_test.shape, Y_test.shape)
In [209]:
model=Sequential()
model.add(LSTM(100,return_sequences=True,input_shape=(30,1)))
model.add(LSTM(100,return_sequences=True))
model.add(LSTM(100))
model.add(Dense(1))
model.compile(loss='mean_squared_error',optimizer='adam')
In [210]:
model.summary()
In [211]:
history = model.fit(X_train, Y_train, epochs=20, batch_size=70, validation_data=(X_test, Y_test),
callbacks=[EarlyStopping(monitor='val_loss', patience=20)], verbose=1, shuffle=False)
In [212]:
train_predict= model.predict(X_train)
test_predict= model.predict(X_test)
#invert predictions
train_predict= scaler.inverse_transform(train_predict)
Y_train=scaler.inverse_transform([Y_train])
test_predict= scaler.inverse_transform(test_predict)
Y_test= scaler.inverse_transform([Y_test])
print('Train Mean Absolute Error:', mean_absolute_error(Y_train[0], train_predict[:,0]))
print('Train Root Mean Squared Error:',np.sqrt(mean_squared_error(Y_train[0], train_predict[:,0])))
print('Test Mean Absolute Error:', mean_absolute_error(Y_test[0], test_predict[:,0]))
print('Test Root Mean Squared Error:',np.sqrt(mean_squared_error(Y_test[0], test_predict[:,0])))
In [213]:
plt.figure(figsize=(8,4))
plt.plot(history.history['loss'], label='Train Loss')
plt.plot(history.history['val_loss'], label='Test Loss')
plt.title('model loss')
plt.ylabel('loss')
plt.xlabel('epochs')
plt.legend(loc='upper right')
plt.show();
In [214]:
aa=[x for x in range(60)]
plt.figure(figsize=(8,4))
plt.plot(aa, Y_test[0][:60], marker='.', label="actual")
plt.plot(aa, test_predict[:,0][:60], 'r', label="prediction")
# plt.tick_params(left=False, labelleft=True) #remove ticks
plt.tight_layout()
sns.despine(top=True)
plt.subplots_adjust(left=0.07)
plt.ylabel('Price', size=15)
plt.xlabel('Time step', size=15)
plt.legend(fontsize=15)
plt.show();
In [215]:
aa=[x for x in range(300)]
plt.figure(figsize=(8,4))
plt.plot(aa, Y_train[0][:300], marker='.', label="actual")
plt.plot(aa, train_predict[:,0][:300], 'r', label="prediction")
# plt.tick_params(left=False, labelleft=True) #remove ticks
plt.tight_layout()
sns.despine(top=True)
plt.subplots_adjust(left=0.07)
plt.ylabel('Price', size=15)
plt.xlabel('Time step', size=15)
plt.legend(fontsize=15)
plt.show();
In [216]:
from sklearn.metrics import r2_score
print('Train R2 Score: ', r2_score(Y_train[0], train_predict[:,0]))
print('Test R2 Score: ', r2_score(Y_test[0], test_predict[:,0]))
In [217]:
len(test)
Out[217]:
In [218]:
x_input=test[70:].reshape(1,-1)
x_input.shape
Out[218]:
In [219]:
temp_input=list(x_input)
temp_input=temp_input[0].tolist()
In [220]:
# demonstrate prediction for next 10 days
from numpy import array
lst_output=[]
n_steps=30
i=0
while(i<30):
if(len(temp_input)>30):
#print(temp_input)
x_input=np.array(temp_input[1:])
print("{} day input {}".format(i,x_input))
x_input=x_input.reshape(1,-1)
x_input = x_input.reshape((1, n_steps, 1))
#print(x_input)
yhat = model.predict(x_input, verbose=0)
print("{} day output {}".format(i,yhat))
temp_input.extend(yhat[0].tolist())
temp_input=temp_input[1:]
#print(temp_input)
lst_output.extend(yhat.tolist())
i=i+1
else:
x_input = x_input.reshape((1, n_steps,1))
yhat = model.predict(x_input, verbose=0)
print(yhat[0])
temp_input.extend(yhat[0].tolist())
print(len(temp_input))
lst_output.extend(yhat.tolist())
i=i+1
# print(lst_output)
In [221]:
day_new=np.arange(1,31)
day_pred=np.arange(31,61)
In [222]:
len(df1)
Out[222]:
In [223]:
plt.plot(day_new,scaler.inverse_transform(df1[2976:]))
plt.plot(day_pred,scaler.inverse_transform(lst_output))
Out[223]:
In [224]:
print(Y_test.shape, Y_train.shape, test_predict.shape)
In [225]:
diesel.index[:train_size-31]
Out[225]:
In [226]:
test_predict.shape
Out[226]:
In [227]:
plt.figure(figsize=(25,5))
plt.plot(df['Diesel'][:train_size], label='train', color='#3f5b9c')
plt.plot(df['Diesel'][train_size:], label='test', color='#9c3f5b')
plt.plot(gasoline.index[-69:], test_predict[:,0], label='prediction', color='#9c3f5b', linestyle='--')
plt.title('Stacked LSTM Diesel Results', size=24)
plt.legend()
plt.show()
In [228]:
stlstm_diesel = test_predict[:,0].copy()
In [229]:
len(stlstm_diesel)
Out[229]:
Jet Fuel Stacked LSTM¶
In [230]:
jet.head()
Out[230]:
In [231]:
df1=jet['Jet Fuel']
df1
Out[231]:
In [232]:
scaler= MinMaxScaler(feature_range=(0,1))
df1=scaler.fit_transform(np.array(df1).reshape(-1,1))
df1 #compressed
Out[232]:
In [233]:
len(df1)
Out[233]:
In [234]:
# train_size = int(len(df1) * 0.80)
# test_size = len(df1) - train_size
train_size=2906
train, test = df1[0:train_size,:], df1[train_size:len(df1),:]
In [235]:
print(train.shape, test.shape)
In [236]:
#convert an array of values into a dataset matrix
def create_dataset(dataset, look_back=1):
#creates window pockkets
X,Y = [], []
for i in range(len(dataset)-look_back-1):
a=dataset[i:(i+look_back),0]
X.append(a)
Y.append(dataset[i + look_back,0])
return np.array(X), np.array(Y)
In [237]:
#reshape into X=t and Y=t+1
look_back = 30
X_train, Y_train = create_dataset(train, look_back)
X_test, Y_test= create_dataset(test, look_back)
In [238]:
print(X_train.shape, Y_train.shape, X_test.shape, Y_test.shape)
In [239]:
# reshape input to be [samples, time steps, features] which is required for LSTM
X_train =X_train.reshape(X_train.shape[0],X_train.shape[1] , 1)
X_test = X_test.reshape(X_test.shape[0],X_test.shape[1] , 1)
In [240]:
print(X_train.shape, Y_train.shape, X_test.shape, Y_test.shape)
In [241]:
model=Sequential()
model.add(LSTM(100,return_sequences=True,input_shape=(30,1)))
model.add(LSTM(100,return_sequences=True))
model.add(LSTM(100))
model.add(Dense(1))
model.compile(loss='mean_squared_error',optimizer='adam')
In [242]:
model.summary()
In [243]:
history = model.fit(X_train, Y_train, epochs=20, batch_size=70, validation_data=(X_test, Y_test),
callbacks=[EarlyStopping(monitor='val_loss', patience=20)], verbose=1, shuffle=False)
In [244]:
train_predict= model.predict(X_train)
test_predict= model.predict(X_test)
#invert predictions
train_predict= scaler.inverse_transform(train_predict)
Y_train=scaler.inverse_transform([Y_train])
test_predict= scaler.inverse_transform(test_predict)
Y_test= scaler.inverse_transform([Y_test])
print('Train Mean Absolute Error:', mean_absolute_error(Y_train[0], train_predict[:,0]))
print('Train Root Mean Squared Error:',np.sqrt(mean_squared_error(Y_train[0], train_predict[:,0])))
print('Test Mean Absolute Error:', mean_absolute_error(Y_test[0], test_predict[:,0]))
print('Test Root Mean Squared Error:',np.sqrt(mean_squared_error(Y_test[0], test_predict[:,0])))
In [245]:
plt.figure(figsize=(8,4))
plt.plot(history.history['loss'], label='Train Loss')
plt.plot(history.history['val_loss'], label='Test Loss')
plt.title('model loss')
plt.ylabel('loss')
plt.xlabel('epochs')
plt.legend(loc='upper right')
plt.show();
In [246]:
aa=[x for x in range(60)]
plt.figure(figsize=(8,4))
plt.plot(aa, Y_test[0][:60], marker='.', label="actual")
plt.plot(aa, test_predict[:,0][:60], 'r', label="prediction")
# plt.tick_params(left=False, labelleft=True) #remove ticks
plt.tight_layout()
sns.despine(top=True)
plt.subplots_adjust(left=0.07)
plt.ylabel('Price', size=15)
plt.xlabel('Time step', size=15)
plt.legend(fontsize=15)
plt.show();
In [247]:
aa=[x for x in range(300)]
plt.figure(figsize=(8,4))
plt.plot(aa, Y_train[0][:300], marker='.', label="actual")
plt.plot(aa, train_predict[:,0][:300], 'r', label="prediction")
# plt.tick_params(left=False, labelleft=True) #remove ticks
plt.tight_layout()
sns.despine(top=True)
plt.subplots_adjust(left=0.07)
plt.ylabel('Price', size=15)
plt.xlabel('Time step', size=15)
plt.legend(fontsize=15)
plt.show();
In [248]:
from sklearn.metrics import r2_score
print('Train R2 Score: ', r2_score(Y_train[0], train_predict[:,0]))
print('Test R2 Score: ', r2_score(Y_test[0], test_predict[:,0]))
In [249]:
len(test)
Out[249]:
In [250]:
x_input=test[70:].reshape(1,-1)
x_input.shape
Out[250]:
In [251]:
temp_input=list(x_input)
temp_input=temp_input[0].tolist()
In [252]:
# demonstrate prediction for next 10 days
from numpy import array
lst_output=[]
n_steps=30
i=0
while(i<30):
if(len(temp_input)>30):
#print(temp_input)
x_input=np.array(temp_input[1:])
print("{} day input {}".format(i,x_input))
x_input=x_input.reshape(1,-1)
x_input = x_input.reshape((1, n_steps, 1))
#print(x_input)
yhat = model.predict(x_input, verbose=0)
print("{} day output {}".format(i,yhat))
temp_input.extend(yhat[0].tolist())
temp_input=temp_input[1:]
#print(temp_input)
lst_output.extend(yhat.tolist())
i=i+1
else:
x_input = x_input.reshape((1, n_steps,1))
yhat = model.predict(x_input, verbose=0)
print(yhat[0])
temp_input.extend(yhat[0].tolist())
print(len(temp_input))
lst_output.extend(yhat.tolist())
i=i+1
# print(lst_output)
In [253]:
day_new=np.arange(1,31)
day_pred=np.arange(31,61)
In [254]:
len(df1)
Out[254]:
In [255]:
plt.plot(day_new,scaler.inverse_transform(df1[2976:]))
plt.plot(day_pred,scaler.inverse_transform(lst_output))
Out[255]:
In [256]:
print(Y_test.shape, Y_train.shape, test_predict.shape)
In [257]:
diesel.index[:train_size-31]
Out[257]:
In [258]:
test_predict.shape
Out[258]:
In [259]:
plt.figure(figsize=(25,5))
plt.plot(df['Jet Fuel'][:train_size], label='train', color='#c24e5d')
plt.plot(df['Jet Fuel'][train_size:], label='test', color='#5dc24e')
plt.plot(gasoline.index[-69:], test_predict[:,0], label='prediction', color='#5dc24e', linestyle='--')
plt.title('Stacked LSTM Jet Results', size=24)
plt.legend()
plt.show()
In [260]:
stlstm_jet = test_predict[:,0].copy()
In [261]:
len(stlstm_jet)
Out[261]:
Part 7: ARIMA¶
Gasoline ARIMA¶
In [262]:
# load required libraries
import sys
import statsmodels as ss
import pandas as pd
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.arima_model import ARIMA
import statsmodels.api as sm
In [263]:
gasoline.head()
Out[263]:
In [264]:
gas = gasoline['Gasoline']
gas.head()
Out[264]:
In [265]:
test_stationarity(gas)
In [266]:
#log transform
gas_log = np.log(gas)
gas_log.head()
Out[266]:
In [267]:
test_stationarity(gas_log)
In [268]:
#transform to stationary
gas_log_diff = gas_log.diff(periods=1)
gas_log_diff
Out[268]:
In [269]:
gas_log_diff.dropna(inplace=True)
In [270]:
test_stationarity(gas_log_diff)
In [271]:
len(gas_log_diff)
Out[271]:
In [272]:
# train_size = int(len(gas_log_diff) * 0.80)
train_size = 2905
In [273]:
train = gas_log_diff[:train_size]
test = gas_log_diff[train_size:]
print("Train shape: " + str(len(train)))
print("Test shape: " + str(len(test)))
In [274]:
print(train.head(3), print(test.head(3)))
In [275]:
from statsmodels.graphics.tsaplots import plot_pacf
title = 'PACF: Log+Diff Gasoline'
lags = 40
plot_pacf(train, title=title, lags=lags, zero=False, method=("ols"));
In [276]:
from statsmodels.graphics.tsaplots import plot_acf
# plot ACF function for PH covid-19 deaths
title = 'ACF: Log+Diff Gasoline'
lags = 40
plot_acf(train, title=title, lags=lags, zero=False);
In [277]:
from statsmodels.tsa.arima.model import ARIMA
# fitting a simple ARIMA(1,1,1) model
model = ARIMA(train, order=(1,1,1))
model_result = model.fit()
model_result.summary()
Out[277]:
In [278]:
from sklearn.metrics import mean_squared_error
# evaluate model MSE and RMSE
start = len(train)
end = len(train) + len(test) - 1
arima111_preds = model_result.predict(start=start, end=end, dynamic=False).rename('ARIMA111 Predictions')
arima111_mse = mean_squared_error(test, arima111_preds)
print("ARIMA(1,1,1) MSE: " + str(arima111_mse))
print("ARIMA(1,1,1) RMSE: " + str(np.sqrt(arima111_mse)))
In [279]:
# use auto_arima to automate ARIMA model selection
from pmdarima.arima import auto_arima
In [280]:
model_auto = auto_arima(
train,
exogenous=None, # don't use exogenous variables
m=1, # don't add seasonal cycles as I'm interested only in an ARIMA model
max_order=None, # don't limit total number of non-seasonal AR and MA components
max_p=15, # use maximum of 15 lags for AR and MA components as implied by the ACF plot of the ARIMA(1,1,1) residuals
max_q=15,
max_d=2, # use max integration of 2 to reach stationarity
maxiter=100,
alpha=0.05, # use 0.05 level of significance
njobs=-1, # fit as many as possible models
trend="ct",
start_P=0, # don't use seasonal components as I'm interested only in an ARIMA model
start_Q=0,
start_D=0,
max_P=0,
max_Q=0,
max_D=0)
In [281]:
model_auto.summary()
Out[281]:
In [282]:
model_auto = ARIMA(train, order=(2,0,1))
model_auto_result = model_auto.fit()
model_auto_result.summary()
Out[282]:
In [283]:
# evaluate model MSE and RMSE
start = len(train)
end = len(train) + len(test) - 1
arima201_preds = model_auto_result.predict(start=start, end=end, dynamic=False).rename('ARIMA201 Predictions')
arima201_mse = mean_squared_error(test, arima201_preds)
print("ARIMA(1,1,1) MSE: " + str(arima201_mse))
print("ARIMA(1,1,1) RMSE: " + str(np.sqrt(arima201_mse)))
In [284]:
arima111_preds.head()
Out[284]:
In [285]:
arima201_preds.head()
Out[285]:
In [286]:
#return the predicted values to reverse diff and reverse log to match the original dataset
arima201_r_diff = inv_diff(gas_log[train_size:], arima201_preds, 1)
arima201_r_diff
Out[286]:
In [287]:
arima201_r_diff_exp = np.exp(arima201_r_diff)
arima201_r_diff_exp
Out[287]:
In [288]:
len(arima201_r_diff_exp)
Out[288]:
In [289]:
plt.figure(figsize=(25,5))
plt.plot(gasoline['Gasoline'][:train_size], label='train', color='#e58a3d')
plt.plot(gasoline['Gasoline'][train_size:], label='test', color='#e53d98')
plt.plot(gasoline.index[-100:], arima201_r_diff_exp, label='arima201', color='#e53d98', linestyle='--')
plt.title('ARIMA(2,0,1) Gasoline Results Inverse Diff+Exp', size=15)
plt.legend()
plt.show()
In [290]:
arima_gas = arima201_r_diff_exp
In [291]:
len(arima_gas)
Out[291]:
In [292]:
print('Test Mean Absolute Error:', mean_absolute_error(test, arima_gas))
print('Test Root Mean Squared Error:',np.sqrt(mean_absolute_error(test, arima_gas)))
print('Test R2 Score: ', r2_score(test, arima_gas))
In [292]:
Diesel ARIMA¶
In [293]:
diesel.head()
Out[293]:
In [294]:
diesel = diesel['Diesel']
diesel.head()
Out[294]:
In [295]:
test_stationarity(diesel)
In [296]:
#log transform
d_log = np.log(diesel)
d_log.head()
Out[296]:
In [297]:
test_stationarity(d_log)
In [298]:
#transform to stationary
d_log_diff = d_log.diff(periods=1)
d_log_diff
Out[298]:
In [299]:
d_log_diff.dropna(inplace=True)
In [300]:
test_stationarity(d_log_diff)
In [301]:
len(d_log_diff)
Out[301]:
In [302]:
# train_size = int(len(d_log_diff) * 0.80)
train_size = 2905
In [303]:
train = d_log_diff[:train_size]
test = d_log_diff[train_size:]
print("Train shape: " + str(len(train)))
print("Test shape: " + str(len(test)))
In [304]:
print(train.head(3), print(test.head(3)))
In [305]:
from statsmodels.graphics.tsaplots import plot_pacf
title = 'PACF: Log+Diff Diesel'
lags = 40
plot_pacf(train, title=title, lags=lags, zero=False, method=("ols"));
In [306]:
from statsmodels.graphics.tsaplots import plot_acf
title = 'ACF: Log+Diff Diesel'
lags = 40
plot_acf(train, title=title, lags=lags, zero=False);
In [307]:
from statsmodels.tsa.arima.model import ARIMA
# fitting a simple ARIMA(1,1,1) model
model = ARIMA(train, order=(1,1,1))
model_result = model.fit()
model_result.summary()
Out[307]:
In [308]:
from sklearn.metrics import mean_squared_error
# evaluate model MSE and RMSE
start = len(train)
end = len(train) + len(test) - 1
arima111_preds = model_result.predict(start=start, end=end, dynamic=False).rename('ARIMA111 Predictions')
arima111_mse = mean_squared_error(test, arima111_preds)
print("ARIMA(1,1,1) MSE: " + str(arima111_mse))
print("ARIMA(1,1,1) RMSE: " + str(np.sqrt(arima111_mse)))
In [309]:
# use auto_arima to automate ARIMA model selection
from pmdarima.arima import auto_arima
In [310]:
model_auto = auto_arima(
train,
exogenous=None, # don't use exogenous variables
m=1, # don't add seasonal cycles as I'm interested only in an ARIMA model
max_order=None, # don't limit total number of non-seasonal AR and MA components
max_p=15, # use maximum of 15 lags for AR and MA components as implied by the ACF plot of the ARIMA(1,1,1) residuals
max_q=15,
max_d=2, # use max integration of 2 to reach stationarity
maxiter=100,
alpha=0.05, # use 0.05 level of significance
njobs=-1, # fit as many as possible models
trend="ct",
start_P=0, # don't use seasonal components as I'm interested only in an ARIMA model
start_Q=0,
start_D=0,
max_P=0,
max_Q=0,
max_D=0)
In [311]:
model_auto.summary()
Out[311]:
In [312]:
model_auto = ARIMA(train, order=(3,0,7))
model_auto_result = model_auto.fit()
model_auto_result.summary()
Out[312]:
In [313]:
# evaluate model MSE and RMSE
start = len(train)
end = len(train) + len(test) - 1
arima307_preds = model_auto_result.predict(start=start, end=end, dynamic=False).rename('ARIMA307 Predictions')
arima307_mse = mean_squared_error(test, arima307_preds)
print("ARIMA(1,1,1) MSE: " + str(arima307_mse))
print("ARIMA(1,1,1) RMSE: " + str(np.sqrt(arima307_mse)))
In [314]:
arima307_preds.head()
Out[314]:
In [315]:
arima111_preds.head()
Out[315]:
In [316]:
#return the predicted values to reverse diff and reverse log to match the original dataset
arima307_r_diff = inv_diff(gas_log[train_size:], arima307_preds, 1)
arima307_r_diff
Out[316]:
In [317]:
arima307_r_diff_exp = np.exp(arima307_r_diff)
arima307_r_diff_exp
Out[317]:
In [318]:
len(arima307_r_diff_exp)
Out[318]:
In [319]:
plt.figure(figsize=(25,5))
plt.plot(df['Diesel'][:train_size], label='train', color='#3f5b9c')
plt.plot(df['Diesel'][train_size:], label='test', color='#9c3f5b')
plt.plot(diesel.index[-100:], arima307_r_diff_exp, label='arima307', color='#9c3f5b', linestyle='--')
plt.title('ARIMA(3,0,7) Diesel Results Inverse Diff+Exp', size=15)
plt.legend()
plt.show()
In [320]:
arima_diesel = arima307_r_diff_exp
len(arima_diesel)
Out[320]:
In [321]:
print('Test Mean Absolute Error:', mean_absolute_error(test, arima_diesel))
print('Test Root Mean Squared Error:',np.sqrt(mean_absolute_error(test, arima_diesel)))
print('Test R2 Score: ', r2_score(test, arima_diesel))
Jet Fuel ARIMA¶
In [322]:
jet.head()
Out[322]:
In [323]:
jet = jet['Jet Fuel']
jet.head()
Out[323]:
In [324]:
test_stationarity(jet)
In [325]:
#log transform
j_log = np.log(jet)
j_log.head()
Out[325]:
In [326]:
test_stationarity(j_log)
In [327]:
#transform to stationary
j_log_diff = j_log.diff(periods=1)
j_log_diff
Out[327]:
In [328]:
j_log_diff.dropna(inplace=True)
In [329]:
test_stationarity(j_log_diff)
In [330]:
len(j_log_diff)
Out[330]:
In [331]:
# train_size = int(len(d_log_diff) * 0.80)
train_size = 2905
In [332]:
train = j_log_diff[:train_size]
test = j_log_diff[train_size:]
print("Train shape: " + str(len(train)))
print("Test shape: " + str(len(test)))
In [333]:
print(train.head(3), print(test.head(3)))
In [334]:
from statsmodels.graphics.tsaplots import plot_pacf
title = 'PACF: Log+Diff Jet Fuel'
lags = 40
plot_pacf(train, title=title, lags=lags, zero=False, method=("ols"));
In [335]:
from statsmodels.graphics.tsaplots import plot_acf
title = 'ACF: Log+Diff Jet Fuel'
lags = 40
plot_acf(train, title=title, lags=lags, zero=False);
In [336]:
from statsmodels.tsa.arima.model import ARIMA
# fitting a simple ARIMA(1,1,1) model
model = ARIMA(train, order=(1,1,1))
model_result = model.fit()
model_result.summary()
Out[336]:
In [337]:
from sklearn.metrics import mean_squared_error
# evaluate model MSE and RMSE
start = len(train)
end = len(train) + len(test) - 1
arima111_preds = model_result.predict(start=start, end=end, dynamic=False).rename('ARIMA111 Predictions')
arima111_mse = mean_squared_error(test, arima111_preds)
print("ARIMA(1,1,1) MSE: " + str(arima111_mse))
print("ARIMA(1,1,1) RMSE: " + str(np.sqrt(arima111_mse)))
In [338]:
# use auto_arima to automate ARIMA model selection
from pmdarima.arima import auto_arima
In [339]:
model_auto = auto_arima(
train,
exogenous=None, # don't use exogenous variables
m=1, # don't add seasonal cycles as I'm interested only in an ARIMA model
max_order=None, # don't limit total number of non-seasonal AR and MA components
max_p=15, # use maximum of 15 lags for AR and MA components as implied by the ACF plot of the ARIMA(1,1,1) residuals
max_q=15,
max_d=2, # use max integration of 2 to reach stationarity
maxiter=100,
alpha=0.05, # use 0.05 level of significance
njobs=-1, # fit as many as possible models
trend="ct",
start_P=0, # don't use seasonal components as I'm interested only in an ARIMA model
start_Q=0,
start_D=0,
max_P=0,
max_Q=0,
max_D=0)
In [340]:
model_auto.summary()
Out[340]:
In [341]:
model_auto = ARIMA(train, order=(0,0,3))
model_auto_result = model_auto.fit()
model_auto_result.summary()
Out[341]:
In [342]:
# evaluate model MSE and RMSE
start = len(train)
end = len(train) + len(test) - 1
arima003_preds = model_auto_result.predict(start=start, end=end, dynamic=False).rename('ARIMA003 Predictions')
arima003_mse = mean_squared_error(test, arima003_preds)
print("ARIMA(0,0,3) MSE: " + str(arima003_mse))
print("ARIMA(0,0,3) RMSE: " + str(np.sqrt(arima003_mse)))
In [343]:
arima003_preds.head()
Out[343]:
In [344]:
arima111_preds.head()
Out[344]:
In [345]:
#return the predicted values to reverse diff and reverse log to match the original dataset
arima003_r_diff = inv_diff(j_log[train_size:], arima003_preds, 1)
arima003_r_diff
Out[345]:
In [346]:
arima003_r_diff_exp = np.exp(arima003_r_diff)
arima003_r_diff_exp
Out[346]:
In [347]:
len(arima003_r_diff_exp)
Out[347]:
In [348]:
plt.figure(figsize=(25,5))
plt.plot(df['Jet Fuel'][:train_size], label='train', color='#c24e5d')
plt.plot(df['Jet Fuel'][train_size:], label='test', color='#5dc24e')
plt.plot(jet.index[-100:], arima003_r_diff_exp, label='arima003', color='#5dc24e', linestyle='--')
plt.title('ARIMA(0,0,3) Jet Fuel Results Inverse Diff+Exp', size=15)
plt.legend()
plt.show()
In [349]:
arima_jet = arima003_r_diff_exp
len(arima_jet)
Out[349]:
In [350]:
print('Test Mean Absolute Error:', mean_absolute_error(test, arima_jet))
print('Test Root Mean Squared Error:',np.sqrt(mean_absolute_error(test, arima_jet)))
print('Test R2 Score: ', r2_score(test, arima_jet))
Part 8: FB Prophet¶
Gasoline Prophet¶
In [351]:
# !pip install fbprophet
In [352]:
import io, os, sys, setuptools, tokenize
from prophet import Prophet
In [353]:
df.reset_index(inplace=True)
df.head()
Out[353]:
In [354]:
# Renaming the column & index for Prophet compatability
gas_prop= df[['Date', 'Gasoline']]
In [355]:
gas_prop.rename(columns = {'Date':'ds', 'Gasoline':'y'}, inplace = True)
gas_prop.head()
Out[355]:
In [356]:
# train_size = int(len(gasoline1) * 0.80)
train_size = 2906
In [357]:
train = gas_prop[:train_size]
test = gas_prop[train_size:]
print("Train shape: " + str(len(train)))
print("Test shape: " + str(len(test)))
In [358]:
test.head()
Out[358]:
In [359]:
# fit model
m = Prophet()
m.fit(train)
Out[359]:
In [360]:
# populate forecast
forecast = m.predict(test)
forecast[['ds', 'yhat', 'yhat_lower', 'yhat_upper']].tail()
Out[360]:
In [361]:
forecast.shape
Out[361]:
In [362]:
# plot forecast
m.plot(forecast);
In [363]:
# plot individual components of forecast: trend, weekly/yearly seasonality,
m.plot_components(forecast);
In [364]:
plt.figure(figsize=(25,5))
plt.plot(df['Date'][:train_size],df['Gasoline'][:train_size], label='train', color='#e58a3d')
plt.plot(df['Date'][train_size:],df['Gasoline'][train_size:], label='test', color='#e53d98')
plt.plot(df['Date'][train_size:], forecast.yhat, label='prediction', color='#e53d98', linestyle='--')
plt.title('Prophet Gasoline Results', size=15)
plt.legend()
plt.show()
In [365]:
gas_prophet = forecast.yhat
len(gas_prophet)
Out[365]:
In [366]:
test
Out[366]:
In [367]:
print('Test Mean Absolute Error:', mean_absolute_error(test.y, gas_prophet))
print('Test Root Mean Squared Error:',np.sqrt(mean_absolute_error(test.y, gas_prophet)))
print('Test R2 Score: ', r2_score(test.y, gas_prophet))
Diesel Prophet¶
In [368]:
# !pip install fbprophet
In [369]:
df.head()
Out[369]:
In [370]:
# Renaming the column & index for Prophet compatability
die_prop= df[['Date', 'Diesel']]
In [371]:
die_prop.rename(columns = {'Date':'ds', 'Diesel':'y'}, inplace = True)
die_prop.head()
Out[371]:
In [372]:
# train_size = int(len(gasoline1) * 0.80)
train_size = 2906
In [373]:
train = die_prop[:train_size]
test = die_prop[train_size:]
print("Train shape: " + str(len(train)))
print("Test shape: " + str(len(test)))
In [374]:
test.head()
Out[374]:
In [375]:
# fit model
m = Prophet()
m.fit(train)
Out[375]:
In [376]:
# populate forecast
forecast = m.predict(test)
forecast[['ds', 'yhat', 'yhat_lower', 'yhat_upper']].tail()
Out[376]:
In [377]:
forecast.shape
Out[377]:
In [378]:
# plot forecast
m.plot(forecast);
In [379]:
# plot individual components of forecast: trend, weekly/yearly seasonality,
m.plot_components(forecast);
In [380]:
plt.figure(figsize=(25,5))
plt.plot(df['Date'][:train_size],df['Diesel'][:train_size], label='train', color='#3f5b9c')
plt.plot(df['Date'][train_size:],df['Diesel'][train_size:], label='test', color='#9c3f5b')
plt.plot(df['Date'][train_size:], forecast.yhat, label='prediction', color='#9c3f5b', linestyle='--')
plt.title('Prophet Diesel Results', size=15)
plt.legend()
plt.show()
In [381]:
diesel_prophet = forecast.yhat
len(diesel_prophet)
Out[381]:
In [382]:
print('Test Mean Absolute Error:', mean_absolute_error(test.y, diesel_prophet))
print('Test Root Mean Squared Error:',np.sqrt(mean_absolute_error(test.y, diesel_prophet)))
print('Test R2 Score: ', r2_score(test.y, diesel_prophet))
Jet Fuel Prophet¶
In [383]:
# !pip install fbprophet
In [384]:
df.head()
Out[384]:
In [385]:
# Renaming the column & index for Prophet compatability
jet_prop= df[['Date', 'Jet Fuel']]
In [386]:
jet_prop.rename(columns = {'Date':'ds', 'Jet Fuel':'y'}, inplace = True)
jet_prop.head()
Out[386]:
In [387]:
# train_size = int(len(gasoline1) * 0.80)
train_size = 2906
In [388]:
train = jet_prop[:train_size]
test = jet_prop[train_size:]
print("Train shape: " + str(len(train)))
print("Test shape: " + str(len(test)))
In [389]:
test.head()
Out[389]:
In [390]:
# fit model
m = Prophet()
m.fit(train)
Out[390]:
In [391]:
# populate forecast
forecast = m.predict(test)
forecast[['ds', 'yhat', 'yhat_lower', 'yhat_upper']].tail()
Out[391]:
In [392]:
forecast.shape
Out[392]:
In [393]:
# plot forecast
m.plot(forecast);
In [394]:
# plot individual components of forecast: trend, weekly/yearly seasonality,
m.plot_components(forecast);
In [395]:
plt.figure(figsize=(25,5))
plt.plot(df['Date'][:train_size],df['Jet Fuel'][:train_size], label='train', color='#c24e5d')
plt.plot(df['Date'][train_size:],df['Jet Fuel'][train_size:], label='test', color='#5dc24e')
plt.plot(df['Date'][train_size:], forecast.yhat, label='prediction', color='#5dc24e', linestyle='--')
plt.title('Prophet Jet Fuel Results', size=15)
plt.legend()
plt.show()
In [396]:
jet_prophet = forecast.yhat
len(jet_prophet)
Out[396]:
In [397]:
print('Test Mean Absolute Error:', mean_absolute_error(test.y, jet_prophet))
print('Test Root Mean Squared Error:',np.sqrt(mean_absolute_error(test.y, jet_prophet)))
print('Test R2 Score: ', r2_score(test.y, jet_prophet))
Part 9: Overall Plot Comparison¶
In [398]:
plt.figure(figsize=(25,10))
# plt.plot(df['Date'][:train_size],df['Gasoline'][:train_size], label='train', color='#e58a3d')
plt.plot(df['Date'][train_size:],df['Gasoline'][train_size:], label='actual', color='#e58a3d', linewidth=3)
plt.plot(df['Date'][-69:], lstm_gas, label='lstm', linestyle='--', linewidth=2)
plt.plot(df['Date'][-69:], stlstm_gas, label='stacked lstm', linestyle='--', linewidth=2)
plt.plot(df['Date'][-100:], arima_gas, label='arima', linestyle='--', linewidth=2)
plt.plot(df['Date'][-100:], gas_prophet, label='fbprophet', linestyle='--', linewidth=2)
plt.title('All Models Gasoline', size=15)
plt.legend()
plt.show()
In [399]:
plt.figure(figsize=(25,10))
# plt.plot(df['Date'][:train_size],df['Diesel'][:train_size], label='train', color='#3f5b9c')
plt.plot(df['Date'][train_size:],df['Diesel'][train_size:], label='test', color='#3f5b9c', linewidth=3)
plt.plot(df['Date'][-69:], lstm_diesel, label='lstm', linestyle='--', linewidth=2)
plt.plot(df['Date'][-69:], stlstm_diesel, label='stacked lstm', linestyle='--', linewidth=2)
plt.plot(df['Date'][-100:], arima_diesel, label='arima', linestyle='--', linewidth=2)
plt.plot(df['Date'][-100:], diesel_prophet, label='fbprophet', linestyle='--', linewidth=2)
plt.title('All Models Diesel', size=15)
plt.legend()
plt.show()
In [400]:
plt.figure(figsize=(25,10))
# plt.plot(df['Date'][:train_size],df['Jet Fuel'][:train_size], label='train', color='#c24e5d')
plt.plot(df['Date'][train_size:],df['Jet Fuel'][train_size:], label='test', color='#c24e5d', linewidth=3)
plt.plot(df['Date'][-69:], lstm_jet, label='lstm', linestyle='--', linewidth=2)
plt.plot(df['Date'][-69:], stlstm_jet, label='stacked lstm', linestyle='--', linewidth=2)
plt.plot(df['Date'][-100:], arima_jet, label='arima', linestyle='--', linewidth=2)
plt.plot(df['Date'][-100:], jet_prophet, label='fbprophet', linestyle='--', linewidth=2)
plt.title('All Models Jet Fuel', size=15)
plt.legend()
plt.show()
